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ABSTRACT 


The  onset  of  Marangoni  convection  in  the  float  zone  of  liquid  silicon  is 
studied  from  a  state  at  rest  in  the  absence  of  gravity.  This  time-dependent 
flow  problem  is  solved  numerically  with  the  aid  of  the  Navier-Stokes 
equations  for  an  axisymmetric  flow  with  nonlinear  free  surface  conditions. 
On  this  free  surface  the  temperature  gradient  is  generated  by  heat  transfer 
and  radiation  from  a  heater,  which  is  symmetrically  located  between  the  two 
walls  of  the  float  zone.  After  a  certain  time,  the  flow  is  asymmetrically 
disturbed  by  moving  the  heater  for  a  short  time  away  from  its  symmetric 
position  and  back.  Three  different  Marangoni  numbers  (based  on  the 
temperature  difference  between  heater  and  melting  point  of  silicon)  are 
computed:  10,400,  30,225,  and  50,050.  The  results  show  that  for  Ma 
—  10,400  the  flow  is  steady  and  stable.  For  the  two  higher  Marangoni 
numbers,  however,  the  disturbed  flows  become  unstable,  and  persistent 
oscillatory  modes  of  0.22  Hz  for  Ma  -  30,225  and  0.27  Hz  for 
Ma  =  50,050  develop.  The  free  surface  itself  is  little  deformed  but 
computation  with  a  flat  surface  confirms  Kazarinoff  and  Wilkowski’s  1989 
restilt  that  the  flat  surface  suppresses  instability.  It  is  conjectured  that  the 
flat  surface  imposes  a  strong  local  flow  symmetry,  which  has  a  damping 
effect,  and  that,  in  the  case  of  a  deformable  free  surface,  instability  and 
maintenance  of  oscillation  occur  when  the  two  rolls  attached  to  the  comers, 
which  are  formed  by  the  wall  and  the  free  surface,  separate  and  become 
vortices  with  extremal  vorticity.  These  two  vortices  interact  rhythmically  by 
alternating  build-up  and  decay. 


ADMINISTRATIVE  INFORMATION 

This  work  was  supported  by  the  NASA-Microgravity  Program  under  Order  No.  C- 
32007-M  of  the  NASA  Lewis  Research  Center. 

INTRODUCTION 

In  outer  space  the  growth  of  silicon  crystals  with  high  purity  can  be  achieved  in  a 
greatly  reduced  gravity  field,  so  that  the  melt  is  not  disturbed  by  fluid  motion  due  to 
buoyancy.  Containerless  devices,  which  are  proposed  to  avoid  chemical  reactions  between 
melt  and  container,  have  interfaces  between  melt  and  an  inert-gas  environment.  However, 
the  temperature  gradient  between  heater  and  solid  silicon,  when  the  heat  is  transmitted 
across  and  along  the  interface,  causes  a  gradient  of  the  surface  tension  that  generates  a  fluid 
flow  in  the  melt.  This  surface-tension  driven  flow  is  called  Marangoni  convection  or 
thermocapillary  convection  and  can  occur  in  steady  and  unsteady  (non-oscillatory  or 
oscillatory)  modes  that  are  separated  by  a  critical  flow  parameter.1-*  The  latter  unsteady 
flow  is  particularly  undesirable  for  crystal  growing  because  an  unsteady  flow  affects  the 
quality  of  silicon  chips  considerably.  Methods  such  as  rotating  the  float  zone,3  MHD 
effects,4  and  counteracting  jets5  have  been  proposed  to  suppress  this  undesired  fluid  motion. 

The  numerical  study  of  Marangoni  convection  in  microgravity  has  been  covered 
extensively  in  the  literature  for  a  number  of  devices  with  various  degrees  of  model 
sophistication.  Review  papers  and  recent  publications  with  background  references  include 
those  by  Ostrach,6  Linde,  Preisser  et  al.,  Schwabe,9  Young  and  Chait,10  and  the  GAMM 
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Workshop,  Roux.11  (The  GAMM  Workshop  dealt  essentially  with  a  “benchmark”  model 
which  is  different  from  this  paper’s  approach). 

This  paper  is  based  on  the  arrangement,  sketched  in  Fig.  1.  Liquid  silicon  is  held 
together  between  two  coaxial  circular-cylindrical  rods  of  solid  silicon  by  the  surface  tension 
of  the  liquid-gas  interface  in  a  microgravity  environment.  The  liquid  bridge  between  the  two 
rods  is  called  the  float  zone.  The  heat  for  the  melting  process  is  provided  by  a  ring  heater 
which  moves  coaxially  over  the  rod.  During  that  time  solid  silicon  is  molten  and  then  freezes 
to  a  monocrystal.  This  process  is  modeled  according  to  Fig.  2  with  the  following 
simplifications:  The  boundaries  between  molten  and  solid  silicon  are  considered  planar  with 
melting  and  freezing  neglected.  On  the  liquid-gas  interface,  the  effect  of  the  gas  on  the  liquid 
shear  stress  is  neglected,  that  is,  the  interface  is  considered  a  “free  surface.”  The  flow 
region  is  assumed  axisymmetric  and,  together  with  the  ring  heater,  fixed.  However,  no 
symmetry  along  the  center  plane  between  the  two  walls  is  enforced  so  that  an  oscillatory  flow 
over  the  whole  float  zone  can  be  studied.  At  the  time  t  —  0  the  liquid  silicon  has  the  melting 
temperature  T  =  Ty,  and  at  this  instant,  the  heater  is  turned  on  to  a  maximum  temperature 
Ts>  Ty.  The  developing  temperature  gradient  along  the  free  surface  creates  the  stress  for 
the  Marangoni  flow  which  in  time  penetrates  into  the  interior  of  the  float  zone.  For  small 
differences  in  the  temperatures  7#  and  Ty  a  steady  state  will  be  reached  after  a  certain  time. 
Beyond  a  critical  difference  Ty  —  Ty  the  flow  is  expected  to  become  unstable,  and  a  time- 
dependent  mode  (non-oscillatory  or  oscillatory)  in  the  direction  parallel  to  the  axis  might 
appear.  It  is  the  purpose  of  this  paper  to  find  these  two  modes  and  to  analyze  the 
corresponding  flow  and  temperature  fields. 

The  source  of  energy,  which  drives  the  Marangoni  convection,  comes  from  the  ring 
heater  and  is  transmitted  through  the  free  surface.  In  mathematical  terms,  a  temperature 
gradient  (perpendicular  to  the  surface)  is  generated  by  heat  radiation  and  heat  transfer  from 
the  ring  heater.  The  boundary  conditions  on  the  free  surface  and  the  geometric 
representation  of  the  adjacent  immediate  flow  region  are  hence  very  sensitive  to  both 
physical  and  numerical  modeling.  It  was  found  by  Kamotani  et  al. u  and  by  Kazarinoff  and 
Wilkowski13  that  the  nonlinear  free  surface  conditions  should  not  be  linearized  (that  is,  a  flat 
surface  should  not  be  assumed)  as  most  other  researchers  have  done  because  surface 
deformation  would  initiate  instability  from  a  critical  Ty  —  Ty  on.  To  check  their  finding,  the 
full  nonlinear  free-surface  conditions  are  retained  in  this  paper.  It  may  be  mentioned  that 
experiments14  indicate  only  a  slight  deviation  of  the  real  surface  from  the  flat  cylindrical  one, 
a  fact  which  might  not,  however,  exclude  a  sensitivity  of  the  deformable  free  surface  toward 
flow  instability  (a  fact  confirmed  in  this  paper). 

The  most  serious  restrictions  of  this  flow  model  are  axisymmetry  and  planar  end  walls 
of  the  rods,  with  melting  and  freezing  processes  neglected.  To  a  lesser  degree  the 
assumption  of  a  fixed  ring  heater  affects  the  study  of  the  real  float  zone.  The  axisymmetry 
imposed  on  the  flow  prevents  a  full  investigation  of  instability  because  the  analysis  of  the 
various  modes  of  instability  is  restricted  to  that  of  axisymmetric  disturbances.  Experiments 
with  earth  gravity  for  a  much  larger  Prandtl  number  than  that  of  liauid  silicon  have  shown 
that  azimuthal  disturbances  are  more  “dangerous”  than  axial  ones.  An  argument  for  the 
importance  of  axisymmetric  modes  was  given  by  Kazarinoff  and  Wilkowski.  The  definitive 
answer  to  the  question  of  whether  axisymmetric  modes  are  crucial  must  be  obtained  from 
future  studies. 

Planar  end  walls  of  the  rods  can  cause  a  singularity  at  the  triple  point  of  the  solid, 
liquid,  and  gaseous  phases.  This  singularity  was  mentioned  by  Zebib  et  al.  In  reality,  as  it 
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is  conjectured  from  Fig.  1,  a  cusp  might  develop  which  prevents  such  a  singularity.  Here,  the 
singularity  cannot  be  avoided  unless  the  assumption  of  a  planar  end  wall  is  abandoned.  Most 
researchers  ignore  this  singularity  because  it  appears  to  have  only  a  local  effect. 

Liquid  silicon  has  a  Prandtl  number  of  0.023.  There  are  only  a  few  papers  which  deal 
with  such  a  low  Prandtl  number.  Non-axisymmetric  oscillations  of  the  float  zone  have  been 
observed  experimentally  for  high  Prandtl-number  flows  by  Kamotani  et  al. 12  who  also  argue 
that  the  deformable  free  surface  is  essential  for  the  onset  of  instability.  Kazarinoff  and 
Wilkowski13  found  numerically  axisymmetric  oscillations  in  a  liquid-silicon  float  zone  and 
Fowlis  and  Roberts18  computed  axisymmetric  oscillations  in  a  rotating  liquid-silicon  float 
zone.  In  a  numerical  study  with  a  general  three-dimensional  flow  code,  Rupp  et  al. 19  found 
azimuthal  periodicity.  However,  their  results  are  restricted  to  a  half-zone,  that  is,  they 
considered  only  the  region  — L  fl  <  z  <  0  with  a  solid  wall  at  z  =  0  in  Fig.  2  and  an 
adiabatic  free  surface.  The  restriction  to  a  half-zone  applies  also  to  recent  instability 
studies.20,21  The  only  study,  the  authors  are  aware  of,  that  was  made  with  almost  the  same 
flow  model  but  with  a  different  numerical  scheme  and  with  a  much  coarser  grid  and  time 
step,  is  recorded  in  a  letter  by  Kazarinoff  and  Wilkowski.13  Unfortunately,  there  is  not 
enough  information  to  make  a  comparison  with  the  present  results,  and  an  investigation  of 
the  discrepancies  is  not  possible. 

THE  INITIAL-BOUNDARY  VALUE  PROBLEM 

The  mathematical  description  of  the  flow  model  is  based  on  the  Navier-Stokes 
equations  for  an  incompressible  Newtonian  fluid.  With  cylindrical  polar  coordinates  r,  <f>,  z 


and  the  corresponding  velocity  components  u,  v,  w,  with  axi symmetry  d/d^>  =  0,  and  v  =  0, 
the  equations  of  motion  and  the  energy  equation  are: 

Ur  +  J+Wt  =  0  ,  (1) 

u,  +  ^-(ru2)r  +  (uw\  =  -^-pr+v[uTT  +  (-“)r  +  u^)  ,  (2) 

r  p  r 

w,  +  j{ruw)r  +  (w2)t  =  -jPt+V  [Wrr  +  jWr  +  w„]  ,  (3) 

T,  +  uTr  +  wTt  =  K[T„  +  jTr  +  Tu]  .  (4) 


Here,  p,  p,  u,  and  K  are  the  pressure,  the  density,  the  kinematic  viscosity,  and  the  thermal 
diffusivity  of  the  liquid  silicon,  respectively.  It  may  be  noticed  that  the  conservation  form  of 
the  nonlinear  terms  in  Eqs.  (2)  and  (3)  is  not  retained  in  Eq.  (4)  for  numerical  reasons. 

The  boundary  conditions  for  the  float  zone  are  (Fig.  2): 

Cold  surfaces: 

Z-±  y,  0  <r  <R:  u  =  w  =pz  =0,  T  =  Tu  .  (5) 
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Axis: 


-  Y  <  2  <  +  J  ,  r  =  0:  «  -  °,  HV  =p,  =  Tr  =  0  .  (6) 

The  free  surface  is  described  by  r  =  h(z,  t),  and  the  boundary  conditions  are10,22: 

ht=  u  -whz  ,  (7) 

p  -  2fjur  +  ti(uz  +  wr)hz  =  -  aw(-jp  +  ~~)  -  oT  Vl  +  /t2  ,  (8) 

*1  K2 

(P  -  2favt)ht  +  +  «v)  =  -  +  °1  Vl  +h\  ,  (9) 

K1  “2 

where  Eqs.  (8)  and  (9)  are  the  stress  components  in  the  r-  and  z-directions.  It  is 


<*  —  —  'I  (J  —  Tu)  .  (11) 

The  constant  coefficient  of  the  surface  tension  is  oy,  7  is  the  temperature  coefficient  of  the 
surface  tension  a,  p  is  the  dynamic  viscosity  with  p  =  pv,  and  1/Rj  and  1/R2  are  the  surface 
curvatures  in  the  planes  z  —  const  and  <t>  =  const ,  respectively. 

For  a  flat  free  surface  with  hz  —  0,  Eq.  (7)  reduces  to  u  —  0,  and  one  obtains  from 
Eqs.  (8)  and  (9)  for  the  shear  stress  at  the  free  surface 

-  0*®  .  (12) 

The  occurrence  of  a  jump  in  wr  at  z  *=  ±L  / 2  can  be  seen  immediately  since  wr  is  equal  to 
zero  along  the  wall  but  nonzero  at  the  triple  point  when  this  point  is  approached  from  the 
free  surface.  wT  can  also  be  interpreted  as  negative  vorticity.  Because  of  the  assumptions  of 
axisymmetry  and  v  =  0,  only  the  ^-component  of  the  vorticity  vector,  defined  by 
=  uz  —  wr,  is  nonzero. 

The  temperature  distribution  of  the  free  surface  is  prescribed  by  the  temperature 
gradient  generated  by  the  heat  conduction  and  radiation  from  the  heater 

,  (13) 


9  *  Tu  +  On  —  Tu)  «p  [• 


(M) 


The  constant  k  is  the  thermal  conductivity,  a  the  heat  transfer  coefficient,  s  the 
Stefan-Boltzmann  constant,  e  the  emissivity,  and  a  is  a  distribution  parameter  of  the  heater 
function  9  which  approximates  the  heater.  Ty  refers  to  the  maximum  temperature  of  the 
heater. 

Equation  (4)  imposes  a  feedback  mechanism  on  the  whole  system  in  the  form  that,  at 
each  time  step,  Eq.  (4)  provides  an  update  of  the  surface  temperature  which  in  turn  affects 
the  boundary  conditions  of  the  flow  field.  The  update  of  the  flow  field  then  modifies  the 
temperature  field,  and  a  new  cycle  starts.  To  improve  the  numerical  stability  the  temperature 
field  is  computed  with  the  known  velocity  field  of  the  previous  time  step.  The  numerical 
error  due  to  this  procedure  is  extremely  small.  Thus,  only  one  cycle  between  the 
temperature  and  velocity  fields  is  necessary. 

As  initial  conditions  it  is  assumed  that  the  fluid  is  at  rest  at  the  melting  temperature 
Ty,  that  p  =  ay,  that  the  free  surface  is  flat  (h  —  R ),  and  that  over  an  extremly  short  time 
span  the  temperature  Ty  of  the  heater  is  gradually  turned  on.  In  order  to  induce  instability 
of  the  float  zone,  which  so  far  is  symmetric  with  respect  to  the  plane  z  —  0,  flow  and 
temperature  fields  can  be  disturbed  at  a  given  instant  by  dislocating  the  ring  heater  for  a 
certain  time  span.  The  initial-boundary  value  problem,  defined  by  Eqs.  (1)  through  (11),  (13) 
and  (14),  can  be  made  dimensionless  in  various  ways.  For  instance,  one  could  scale  the 
velocity  either  by  K  /R  or  by  'yAT  jp.  In  the  latter  case  the  following  characteristic  quantities 
would  appear:  length  R,  velocity  'yAt  /p,  time  Rp/'yAT,  temperature  AT,  and  pressure 
(/yAT)2  /pv.  In  addition  to  the  Prandtl  and  Marangoni  numbers  (or  Reynolds  number 
Re  —  Ma /Pr),  the  Weber  number  'yAT fay  (which  is  here  equal  to  the  Capillary  number), 
the  Nusselt  number  aR/k,  and  the  radiation  number  seRAT3/k  will  occur.  Since  no 
simplifications  (truncations)  of  the  original  boundary-value  problem  have  been  made,  the 
dimensional  form  as  well  as  the  dimensionless  models  are  equivalent.  As  in  Refs.  13  and  16, 
because  of  the  importance  of  silicon,  the  dimensional  form  is  retained.  The  following  data 
from  Refs.  3  and  10  are  used: 

Table  1.  List  of  constants  for  silicon* 


p  =  2.5  g/cm3 
p  —  0.0088  g/cm  s 
v  =  0.0035  cmz/s 
k  —  0.32-107  erg/cm  s  °C 
K  =  0.15  cm2/s 
a  =  0.64-105  erg/s  cm 2  °C 
ay  =  720  dyn  /cm 
7=  0.43  dyn /cm  °C 
Ty  =  1685  K 
c  =  0.3 
Pr  *  0.023 

Ma  —  325  RAT,  AT  =  Tu  -  Tu 

The  Stefan-Boltzmann  constant  is  s  =  5.668-10-5  erg  /cm1  s  K4,  and  the  Prandtl  number  Pr 
and  Marangoni  number  Ma  are  defined  by 

>  Ma- •  (15) 
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Another  Marangoni  number  can  be  defined  that  is  based  on  the  temperature  difference 
between  the  highest  surface  temperature  Tmax  and  Tm .  This  Marangoni  number,  designated 
by  Mo*  =  R~i(Tmax  —  TMyfiK,  is  not  an  independent  parameter  but  part  of  the  solution. 

OUTLINE  OF  THE  NUMERICAL  PROCEDURE 

The  numerical  solution  of  the  initial-boundary  value  problem,  defined  by  Eqs.  (1) 
through  (11),  (13)  and  (14),  is  carried  out  with  the  aid  of  a  finite-difference  technique, 
boundary-fitted  coordinates,  and  artificial  compressibility.  Details  are  given  in  Ohring  and 
Lugt,23  with  applications  to  different  flow  problems.24,25  Thus,  only  an  outline  is  given  in  this 
paper.  Hsieh  and  Pline,26  in  a  numerical  scheme  somewhat  similar  to  the  present  one,  also 
employ  boundary-fitted  coordinates  and  artificial  compressibility.  In  addition,  they  use  a 
third-order  upwind  scheme  for  the  convective  terms,  whereas  in  this  paper  central 
differencing  without  added  smoothing  terms  are  used  for  the  convective  terms.  Reviews27,28 
on  numerical  methods  for  free  boundaries  also  describe  other  possible  methods. 

For  the  numerical  integration  of  this  initial-boundary  value  problem  it  is  convenient  to 
make  a  boundary-fitted  coordinate  transformation.  In  the  figures  that  follow,  the  origin  of 
the  z-coordinate  is  shifted  to  — L  /2,  that  is,  the  new  coordinate  is 
z'  —  z  +Lfl=  z  +0.5  cm.  Figure  3a  is  a  schematic  drawing  of  the  way  in  which  the 
physical  plane  (z\  r)  is  mapped  onto  the  computational  domain  ((,  rj).  Only  the  coordinate 
lines  which  form  the  boundaries  of  the  two  regions  are  drawn.  The  coordinate  lines  in 
physical  space  are  mapped  onto  a  uniformly  spaced  Cartesian  mesh  with  a  unit  mesh  spacing 
in  each  coordinate  direction. 

As  the  flow  field  evolves  in  time,  the  grid  in  physical  space  will  move  and  its 
coordinate  lines  will  be  attracted  to  regions  of  high  flow  gradients  through  the  use  of  an 
adaptive-grid  technique.  However,  the  Cartesian  grid  in  computational  space  always  remains 
fixed  and  uniform.  This  is  the  major  advantage  of  using  a  mapping.  The  physical  region  is 
mapped  onto  a  computational  space  with  a  Cartesian  grid  consisting  of  161  points  in  the 
(—coordinate  direction  and  99  points  in  the  rj— coordinate  direction.  Total  mass  conservation 
at  each  time  step  is  enforced  by  distributing  any  deficient  or  surplus  amount  evenly  along  the 
free  surface.  (It  may  be  mentioned  that  the  number  of  grid  points  used  in  this  paper  is 
almost  two  orders  of  magnitude  greater  than  that  of  Ref.  13,  a  fact  of  particular  importance 
for  the  accuracy  of  the  solution  at  the  free  surface). 

The  curvilinear  coordinates  ((,  rj)  are  obtained  as  solutions  of  the  two  elliptic  partial 
differential  equations  with  the  physical  space  coordinates  (z,  r)  as  independent  variables 


{„+{„=  <{?+{?)**«.  1) , 

(16) 

1 kr  +  >fc  -  l)  • 

(17) 

P*  and  QT  are  control  functions.23  Since  all  calculations  are  to  be  done  in  the  rectangular 
computational  domain,  these  two  elliptic  partial  differential  equations  are  transformed  by 
interchanging  the  dependent  and  independent  variables.  As  a  result,  the  physical  space 
coordinates  (z,  r)  are  solved  in  terms  of  the  computational  space  coordinates  ((,  rj)  at  each 
time  step.29 

The  continuity  equation  (1)  is  replaced  by  an  equation  with  pseudo-compressibility  for 
numerically  conserving  mass  at  each  physical  time  step: 
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(18) 


p  t  +  ur  -f  ~  +  iVj  =  0  . 

t  is  the  pseudo-time.  Mass  conservation  was  excellent  for  all  flow  cases  computed: 
Typically,  the  average  value  of  |  V-T?|  for  the  largest  Marangoni  number  was  between  lO-4 
and  10“3  s"1,  except  when  the  flow  was  disturbed  by  the  moving  heater.  Then,  the  average 
value  of  |  V-T*]  was  between  10 ~3  and  10"2  j-1.  These  data  may  be  compared  with  the 
maximum  value  of  vorticity  of  1350  s-1. 

Altogether  six  partial  differential  equations  for  u,  w,  T,  r,  z,  and  p  must  be  solved 
with  the  proper  boundary  conditions.  The  finite-difference  technique  for  solving  these 
equations  is  briefly  described  in  the  following  way:  All  spatial  derivatives,  including  one¬ 
sided  derivatives  at  the  boundaries,  are  replaced  by  finite-difference  operators  of  second 
order  in  the  computational  space.  The  time-differencing  procedure  is  implicit  with  a  certain 
number  of  pseudo-time  steps  for  each  physical  time  step. 

A  “four-color”  scheme  (Fig.  3b)  is  used  in  the  interior  of  the  computational  space. 
The  use  of  such  a  scheme,  which  can  be  vectorized,  resulted  in  an  order  of  magnitude 
increase  in  computer  speed  on  the  Cray-XMP  2/16  and  Cray-YMP  8/64  on  which  the 
computations  were  performed.  The  “four-color”  scheme  consists  of  obtaining  updates  for  u 
and  w  at  all  the  o  points  simultaneously,  then  at  all  the  |  points,  the  x  points,  and  the  A 
points,  in  that  order.  The  latest  available  updates  are  used  in  this  process.  The  “four-color” 
scheme  is  also  applied  to  p  and  to  T  independently. 

The  computational  cycle  for  one  complete  pseudo-time  step  iteration  consists  of  (a) 
applying  the  “four-color”  scheme  to  compute  updates  for  r  and  z  followed  by  obtaining  the 
latest  updates  for  r  at  successive  points  along  the  left  and  right  boundaries;  (b)  obtaining 
updates  for  p  at  successive  points  along  the  boundaries  from  the  boundary  conditions  for  p 
followed  by  applying  the  “four-color”  scheme  to  compute  updates  for  p  in  the  interior;  and 
(c)  obtaining  updates  for  u  and  w  at  successive  points  along  the  boundaries  from  their 
boundary  conditions,  and  then  obtaining  updates  for  h  at  successive  points  along  the  free 
surface  from  the  kinematic  condition  (7),  followed  by  applying  the  “four-color”  scheme  to 
compute  updates  for  u  and  w  in  the  interior. 

At  the  completion  of  this  computational  cycle,  after  the  latest  updates  for  r,  z,  u,  and 
w  satisfy  certain  convergence  criteria  at  all  points,  these  updates  are  tire  solution  at  the  new 
time  level  n  + 1.  If  the  convergence  criteria  are  not  met,  cycle  (a)  through  (c)  is  repeated 
until  they  are  met.  A  new  time  step  is  considered  computed  when  the  corresponding 
components  of  the  velocity  fields  of  two  consecutive  iterations  are  within  1%  of  each  other. 
The  accuracy  of  a  very  similar  numerical  scheme  was  checked  with  fine  grids.23  It  may  be 
mentioned  that  the  singularities  at  r  =  R,  z  —  ±  L  p.  diminish  the  efficiency  of  the 
numerical  scheme  considerably,  but  their  range  of  influence  on  the  whole  flow  field  is  only 
local. 


Before  the  iterative  process  just  described  is  used  for  obtaining  r,  z,  u,  w  and  p  at 
time  level  n  + 1,  the  temperature  T  at  the  new  time  level  n  + 1  is  obtained  by  the  same 
iterative  method  with  the  velocities  u  and  w  from  time  level  n.  Computations  of  the  stream 
function  are  not  presented  in  this  paper.  Therefore,  the  streamline  pictures  show  only 
selected  streamlines  obtained  numerically  from  the  velocity  fields.  These  selected 
streamlines  do  not  represent  equally  spaced  incremental  values  of  the  mass  flux. 
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RESULTS 


Three  numerical  calculations  were  performed  for  Ma  —  10,400,  30,225,  and  50,050 
with  the  constants  of  Table  1,  with  R  =  L  =  1  cm  and  with  a  =  10-4.  The  flow  was 
disturbed  for  all  three  cases  at  t  =  1.52  s  by  moving  the  heater  from  z  =  0  to  z  =  0.25  cm 
at  t  =  1.62  s  and  back  at  time  t  =  1.72  s. 

Before  the  results  are  given,  the  following  definitions  and  remarks  are  useful  with 
regard  to  the  graphical  presentation  of  the  data  and  their  interpretation.  An  axisymmetric 
vortex  ring  is  defined  as  a  fluid  motion  with  nested  closed  streamlines  in  the  meridional  plane 
(r,  z)  with  reference  to  a  frame  fixed  to  the  center  of  the  streamlines.30  In  a  steady  flow,  the 
refeience  frame  is  the  coordinate  system  (r,  z)  itself.  For  the  vorticity  field  in  such  a  steady 
flow  with  closed  streamlines,  Prandtl’s  restriction31  holds  that  the  integral  over  the  viscous 
terms  along  closed  streamlines  vanishes 

v  curl'd  •  dT=  0  ,  (19) 

with  di T  the  length  element  along  the  streamlines.  Prandtl  concluded  that  for  almost 
frictionless  motion,  that  is  for  an  almost  inviscid  fluid  motion  with  arbitrarily  small  v  but 
i/#0,  the  vorticity  of  an  axisymmetric  flow  in  the  meridional  plane  inside  a  closed 
streamline  must  be 


Uj  =  const  •  r  .  (20) 

This  result  follows  from  Helmholtz’s  conservation  law  d(u+/r)/dt  =  0  for  an  axisymmetric 
inviscid  fluid  flow  with  v  —  0  that  yields  for  a  steady-state  flow 

f  =/(*)■  (21) 

rp  is  the  stream  function.  Eq.  (20)  is  a  special  case  of  a  maximum  principle  that  holds  for 
any  steady  axisymmetric  motion  of  a  viscous  fluid:  u+/r  can  have  an  extremum  only  at  the 
boundary.32  A  vortex  ring  (characterized  above  by  nested  closed  streamlines)  in  a  steady 
flow  is  thus  “attached”  in  the  sense  that  its  a )+/r  -  field  has  no  extremum  inside  the  flow 
region  and  that  lines  of  constant  ui^/r  may  at  the  most  form  “tongues.”  However,  this 
vorticity  field  without  extremum  still  must  obey  restriction  (19). 

Vortex  shedding  is  an  intrinsically  time-dependent  process  when  the  tongues  of  u+/r- 
lines  pinch  off  and  form  closed  lines.  This  means  that  an  extremum  of  oj+/r  has  developed 
and  that  the  vortex  ring  has  “detached.”  Because  of  the  importance  of  the  u^/r  -  field, 
vorticity  is  presented  in  the  text  that  follows  in  the  form  of  lines  of  constant  o^/r.  Some 
examples  of  equivorticity  lines,  that  is  of  lines  of  constant  are  also  given  for  comparison. 
These  lines  are  significant  for  the  study  of  the  flux  of  vorticity,  defined  by  vdu+/dn,  with  n  the 
coordinate  normal  to  the  lines.  That  dut^/dn  is  relevant  and  not  d{ui^/r)/dn  can  be  seen 
immediately  in  the  example  of  the  Poiseuille  flow  in  a  pipe,  for  which  d(ui+/r)/dn  is  zero  at 
the  wall  but  not  do o+fdn.  However,  for  all  practical  purposes  the  flux  inside  the  fluid  can  be 
discussed  also  with  lines  of  constant  o^/r. 

It  may  be  mentioned  that  the  center  of  nested  closed  streamlines,  the  location  of 
extremum  of  u^/r,  and  the  location  of  extremum  of  do  not  coincide  in  general.  The 
results  will  show  that  the  center  of  closed  streamlines  is  closer  to  the  location  of  extremum 
of  u+/r  than  to  the  location  of  extremum  of 
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The  Case  Ma  =  10,400 


The  Marangoni  number  Ma  =  10,400  corresponds  to  a  temperature  difference 
Th  -  Tm  of  32°C.  Figure  4  shows  lines  of  constant  u^/r  at  six  different  times  (solid  lines 
represent  negative  and  dashed  lines  positive  contours)  prior  to  and  after  the  asymmetric 
disturbance  of  the  flow.  After  the  heater  is  turned  on  at  /  =  0,  immediately  strong  layers  of 
vorticity  along  the  free  surface  are  generated  by  the  surface  shear  stress  that  increase  in 
strength  toward  the  corners  (which  are  formed  by  the  walls  and  the  free  surface). 
Associated  with  the  heater-generated  surface  shear  stress  is  the  occurrence  of  surface 
velocity  which  increases  toward  the  corners,  builds  up  in  time,  and  reaches  values  up  to 
2.93  cm  Js.  This  process  is  called  slamming.  The  surface  vorticity  is  forced  away  from  the 
walls  and  deflected  into  the  interior  by  forming  tongues.  This  curling  deflection  of  surface 
vorticity  due  to  slamming  creates  wall  vorticity  of  opposite  sign  that  also  forms  tongues 
beneath.  The  flow  situation  is  different  from  ordinary  cavity  flows  in  which  the  sign  of  the 
vorticity  does  not  change  at  the  corner,  and  a  boundary  layer  along  the  entire  solid  wall 
develops.33  In  Fig.  4  no  vortex  shedding,  that  is,  no  extremum  of  w^/r,  is  observed.  At 
t  —  10.5  s  the  flow  reaches  an  almost  steady  state.  The  enlarged  core  of  the  right  roll  is 
displayed  in  Fig.  5  which  shows  lines  of  constant  u^/r  indicating  a  smooth  transition  from 
high  corner  values  to  zero  at  the  center  point  r  =  1,  z'  =  0.5  cm.  Thus,  the  maximum 
principle  for  cj +/r  is  obeyed;  the  roll  is  “attached.”  For  comparison.  Fig.  6  demonstrates 
the  occurrence  of  extrema  of  and  hence  the  unsuitability  of  the  -  field  to  distinguish 
between  steady  state  and  shed  vortices.  In  Fig.  7  the  surface  vorticity  (u^),  is  displayed 
which  is  equal  to  the  negative  shear  stress  at  the  surface.  Since  this  free  surface  is  barely 
deformed  (that  agrees  with  Kazarinoff  and  Wilkowski’s  findings13),  the  surface  vorticity  can 
be  checked,  for  the  computed  Tt  given,  with  Eq.  (12)  for  a  flat  surface.  The  agreement  is  so 
good  that  the  two  curves,  shown  in  Fig.  7  for  the  times  1.00  s,  1.72  s,  and  10.5  s ,  are 
indistinguishable,  except  at  the  extrema  near  the  wall.  Surface  vorticity  increases  with  time 
near  the  walls.  The  tongues  formed  by  the  surface  and  wall  vorticity  prevent  the  two  toroidal 
convection  rolls  to  occupy  the  whole  float  zone  and  confine  the  main  fluid  motion  to  the 
surface  region.  The  streamlines  in  Fig.  8  show  these  rolls  with  the  flow  direction  from 
warmer  to  cooler  parts  of  the  free  surface.  The  center  of  the  rolls  moves  slightly  toward  the 
symmetry  line  z'  =  0.5  cm  but  stays  close  to  the  free  surface.  At  almost  steady  state 
t  =  10.5  s,  two  additional  rolls  beneath  the  main  ones  occur  (they  have  been  observed 
already  at  /  =  7.5  s  but  not  shown  in  Fig.  8). 

The  isotherms  in  Fig.  9  reveal  that  the  temperature  gradient  at  the  cold  walls  increases 
with  time  but  is  unevenly  distributed  along  the  walls.  This  means  that  the  heat  transfer  from 
the  liquid  to  the  solid  phase  increases  with  time  but  much  more  near  the  free  surface.  This 
will  be  verified  quantitatively  further  below.  The  temperature  at  the  free  surface  itself,  that  is 
the  difference  (Tt  —  T^)°C,  is  plotted  in  Fig.  10a  for  times  before  the  disturbance  is 
imposed.  The  curves  end  at  the  cold  wall  at  a  nonzero  angle.  The  corner  points 
(r  —  1,  z'  —  0  and  1  cm)  are  thus  singular.  The  values  of  the  Marangoni  number  Ma*  are 
230  at  /  —  0.50  s,  268  at  /  =  1.0  s,  and  278  at  /  =  1.52  s. 

Despite  the  strong  disturbance  imposed  at  /  =  1.52  s,  there  is  no  tendency  toward 
instability.  The  flow  field  returns  quickly  to  an  almost  steady  state  as  demonstrated  in  Fig.  11 
where  the  surface  vorticity  at  z'  =  0.5  s  is  plotted  against  time.  A  tiny  damped  oscillation  is 
visible  that  is  negligible  if  one  considers  that  at  t  —  4.0  s  the  amplitude  is  less  than  2%  of  the 
maximum  surface  vorticity.  The  lines  of  constant  u+fr  in  Fig.  4  at  /  =  10.5  s  display 
symmetric  tongues  except  for  the  lowest  level  of  vorticity  of  the  order  one  s~x.  This 
asymmetry  appears  only  at  a  level  of  less  than  2%  of  the  maximum  surface  vorticity  which 
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reaches  a  value  of  |  (w^)|  s  =  140  s  1  (Fig.  7).  The  surface  velocity  \vs  for  t  =  10.5  s  is  given 
in  Fig.  12  and  reveals  complete  symmetry. 

The  temperature  fields  in  Fig.  9  change  barely  after  the  disturbance  except  for  a  thin 
layer  along  the  free  surface  as  shown  in  Fig.  10b.  From  t  =  7.5  s  on  the  surface  temperature 
is  symmetric  around  z'  =  0.5  cm  and  the  two  curves  for  t  =  7.5  s  and  t  =  10.5  5  coincide 
completely  with  a  maximum  surface  temperature  difference  of  T,  -  TM  =  0.80  which 
corresponds  to  a  value  of  Ma*  =  260°C.  This  coincidence  confirms  that  the  fluid  motion  is 
at  an  almost  steady  state.  The  result  does  not  agree  with  Kazarinoff  and  Wilkowki’s  finding13 
that  instability  should  occur  at  a  Marangoni  number  of  Ma  =  1202.  In  this  context  it  may  be 
mentioned  that  symmetric  initial  solutions  should  remain  symmetric  in  time  with  an  accurate 
numerical  scheme  and  should  only  develop  asymmetrical  patterns  after  an  asymmetric 
disturbance  is  imposed. 

The  Case  Ma  =  50,050 

The  second  case  with  Ma  =  50,050  corresponds  to  a  temperature  difference  of  154 °C. 
As  in  the  previous  case,  the  heater  was  turned  on  at  t  =  0  in  a  symmetric  position  and  then, 
at  /  =  1.52,  the  heater  was  moved  out  of  this  symmetric  position  and  back  to  create  an 
asymmetric  disturbance  of  the  flow  field.  Lines  of  constant  c^/r,  surface  vorticity,  and 
streamlines  are  shown  in  Figs.  13  through  16,  and  isotherms  and  surface  temperature  in  Figs. 
18  and  19  for  the  period  from  t  =  1.50  s  to  4.00  s. 

As  expected,  the  slamming  of  the  flow  near  the  surface  against  the  cold  walls  is  larger, 
compared  to  the  flow  of  the  lower  Marangoni  number.  The  core  of  the  rolls  is  wider  (Fig. 
13),  and  curved  tongues  of  vorticity  from  the  free  surface  fill  out  the  core,  without  forming 
closed  <jj+/r  -  lines  before  the  disturbance  (Fig,  14,  t  —  1.5  s ).  Fig.  15  shows  the  magnitude 
of  the  surface  vorticity  which  reaches  a  value  of  1200  s-1  at  t  —  2.5  s.  Maximum  surface 
velocities  of  about  10  cm /s  occur.  The  two  additional  weak  toroidal  rolls  beneath  the  main 
rolls  develop  earlier  than  for  Ma  —  10,400  (Fig.  16). 

The  disturbance  imposed  causes  now  an  asymmetric  pair  of  rolls  and  permanent 
oscillation.  The  rolls  have  extrema  of  u^/r  as  seen  in  Fig.  14,  t  =  2.5  s,  indicating  detached 
rolls.  Compared  to  attached  vortices  characterized  by  vorticity  tongues  only,  detached 
vortices  have  a  larger  moment  of  inertia  due  to  higher  core  rotation.  The  way  detached  rolls 
increase  and  then  maintain  the  oscillation  of  the  float  zone  for  Ma  —  50,050  will  be 
explained  in  detail  below.  It  suffices  here  to  demonstrate  the  undamped  oscillation.  In  Fig. 
17  the  surface  vorticity  at  z'  —  0.5  cm  and  z'  =  0.0125  cm  is  plotted  against  time.  The 
amplitude  increases  and  approaches  a  constant  value.  The  frequency  of  the  oscillations  is 
0.27  Hz  and  is  constant  over  the  time  span  computed. 

Isotherms  and  surface  temperature  are  plotted  in  Figs.  18  and  19.  Major  changes  in 
the  temperature  field  occur  only  in  a  layer  adjacent  to  the  free  surface.  The  surface  itself 
deforms  on  a  “microscopic”  scale  which  is  of  the  order  of  10-4  cm,  in  fair  agreement  with 
Kazarinoff  and  Wilkowski.13 

The  last  two  cycles  from  t  —  25  s  to  /  =  34  s  have  almost  constant  amplitude,  and 
their  behavior  is  recorded  in  Figs.  20  through  29.  Lines  of  constant  u^/r  in  Figs.  20  and  21 
and  equivorticity  lines  in  Fig.  22  reveal  the  oscillatory  behavior  of  the  flow  field  in  the  form 
of  rolls  of  pulsating  size.  The  main  activity  is  restricted  to  the  upper  half  of  the  float  zone 
near  the  free  surface.  The  corresponding  surface  vorticity,  surface  velocity,  streamlines, 
isotherms,  and  surface  temperature  are  plotted  respectively  in  Figs.  23  through  27. 
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A  close-up  of  the  core  region  of  the  rolls  in  Fig.  21  is  being  used  now  to  describe  the 
cyclic  change  of  the  float  zone.  During  the  transient  phase,  each  of  the  two  rolls  has 
established  a  core  of  vorticity  with  extreme  values  (visible  as  nested  closed  ^/r-lines  except 
at  certain  transition  phases).  This  local  extremum  assumes  maxima  and  minima  during  a 
cycle  representing  periods  of  accelerating  and  decelerating  core  rotation.  The  alternating 
build-up  and  decay  of  the  rolls  is  accompanied  by  vorticity  flux,  changing  its  direction  and 
visible  in  the  form  of  retreating  and  advancing  vorticity  tongues  and  closed  lines  as  well.  The 
direction  of  the  flux  in  the  right  core  is  indicated  by  arrows  in  Fig.  21.  A  comparison  with 
equivorticity  lines  in  Fig.  22  again  demonstrates  that  u^/r  -  lines  are  better  suited  for 
displaying  the  flow  characteristics  than  u -  lines.  Notice  the  occurrence  of  two  extrema  in 
Fig.  22  at  t  =  27,  28,  and  30  s. 

The  processes  involved  are  described  now  in  detail.  To  begin  with,  the  flow  field  at 
/  =  28  s  in  Fig.  21  is  chosen  when  the  largest  displacement  of  the  flow  (that  is,  the  largest 
amplitude)  from  the  “symmetric”  flow  configuration  occurs.  The  right  roll  has  almost 
attained  its  smallest  size  associated  with  a  maximum  of  |  u>+/r  |  «  |  —  47  J  cm-1s-1  at  the 
core’s  center.  The  roll  on  the  left  side,  in  contrast,  has  acquired  its  largest  extension  and  is 
approaching  its  minimum  of  positive  u^/r.  The  surface  temperature  has  its  highest  value  in 
this  left  half  of  the  float  zone  (Fig.  27)  and  provides  maximum  sunt.,  vorticity  and  slamming 
(Figs.  23  and  24).  Surface  vorticity  is  now  pouring  into  the  center  and  is  filling  up  the  core 
(Fig.  21,  /  —  29  s).  This  means  that  the  rotation  of  the  left  roll  increases  and  the  roll 
contracts.  This  process  takes  place  in  concert  with  the  right  roll  which  expands  through 
decay,  that  is,  through  diffusion  of  center  vorticity.  An  instant  will  be  reached  (/  —  30  s),  at 
which  the  left  roll  has  shrunk  to  its  minimum;  the  center  vorticity  of  the  left  roll  has  leveled 
and  closed  lines  vanish  (Fig.  14,  t  —  4.0  s  shows  also  such  an  instant.  At  the  free  surface, 
meanwhile,  the  location  of  zero  vorticity  shifts  to  the  right  because  of  the  strengthening  of 
surface  vorticity  on  the  left  side.  This,  in  turn,  causes  the  shift  of  maximum  surface 
temperature  to  the  right  since  zero  surface  vorticity  and  maximum  surface  temperature  are 
tied  together  (Fig.  27).  Thus,  despite  the  symmetric  position  of  the  ring  heater,  the  location 
of  maximum  surface  temperature  depends  on  the  fluctuating  flow  field.  After  t  =  30  s  the 
process  repeats  itself  with  exchanged  roles. 

Based  on  this  description  of  the  pulsating  core  region,  the  following  explanation  is 
offered.  About  t  —  28  s  (Fig.  21),  the  right  roll  shrinks  to  a  concentrated  vortex.  The  thin 
feeding  layer  of  vorticity,  beginning  at  the  corner  and  surrounding  the  core,  strengthens 
through  contraction.  At  die  same  time  the  left  roll  expands;  the  surrounding  feeding  layer  is 
stretched  and  weakened.  The  combined  action  of  the  two  vorticity  layers  at  the  surface 
causes  the  location  of  u,  —  0,  that  is  to  move  to  the  left.  The  maximum 

displacement  is  reached  at  about  t  —  28  s  (Fig.  29),  when  the  build-up  and  emptying  of  the 
two  rolls,  respectively,  have  ceased.  Then,  the  process  is  reversed  (Fig.  21,  t  =  29  s).  The 
restoring  mechanism,  which  tries  to  re-establish  the  symmetric  flow  situation  (zero 
amplitude),  consists  now  of  contracting  the  left  roll  and  of  stretching  the  right  roll.  This 
restoring  mechanism  is  supported  by  a  maximum  input  of  surface  vorticity  at  the  left  corner. 
At  /  *=  29  s,  zero  amplitude  occurs  with  an  overshoot  that  is  visible  in  the  non-symmetry  of 
the  flow  (Fig.  21)  and  caused  by  the  moment  of  inertia  of  the  rolls.  The  analogy  to  a 
swinging  pendulum  is  at  hand. 

It  becomes  clear  now  what  happens  during  the  transient  phase  after  the  heater  has 
briefly  been  moved  to  the  right  and  back.  At  /  =  1.72  s,  the  right  roll  is  about  to  form 
closed  vorticity  lines  with  the  largest  rotation  next  to  the  the  center,  while  the  left  roll  has 
still  the  feature  of  an  attached  vortex  (not  shown  in  Fig.  14).  At  t  =  2.5  s  in  Fig.  14,  the 
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right  roll  has  clearly  established  itself  as  a  detached  vortex  with  closed  loops  and  with 
vorticity  pouring  into  the  center.  The  left  roll  is  weaker  but  occupies  yet  about  the  same 
space  as  the  right  roll.  The  stronger  right  roll  has  pushed  (Ts)max  a  little  into  the  left  half.  At 
t  —  4.0  s  the  reverse  situation  is  seen.  With  each  cycle,  a  larger  body  of  fluid  gets  involved 
in  the  oscillation.  The  expanding  roll  requires  more  space  and  squeezes  the  other  roll  which 
in  turn  gets  more  concentrated.  This  process  is  supported  by  the  simultaneous  larger  shift  of 
(r,)BU  away  from  the  symmetry  line  z  =  0.  The  amplitude  of  the  oscillation  increases  until 
the  vorticity  gradients  in  and  around  the  cores  become  so  large  that  diffusion  takes  over.  It 
is  remarkable  that  the  frequency  of  the  oscillation  during  the  transient  period  remains 
constant. 

The  growth  and  maintenance  of  oscillation  for  Ma  —  50,050  is  in  contrast  to  the  flow 
situation  for  Ma  =  10,400.  Here,  the  asymmetrically  disturbed  flow  near  the  surface  is 
immediately  damped  by  diffusion  and  the  core  regions  of  the  roll  barely  affected  (Fie  4). 
There  is  a  certain  analogy  of  the  onset  of  instability  and  the  maintenance  of  the  osc  g 

float  zone  with  the  symmetry  breaking  of  a  parallel  flow  behind  a  circular  cylinder  e 

steady  symmetric  flow  with  attached  vortices  changes  at  a  critical  Reynolds  number  io  a 
periodic  vortex  street  with  detached  vortices.  It  may  be  mentioned  that  at  the  same 
Reynolds  number  beyond  the  critical  one  the  two  (unstable)  attached  vortices  of  the  steady- 
state  configuration  are  weaker  than  the  shed  vortices  of  the  (stable)  unsteady  configuration. 
The  latter  one  is  apparently  more  energy-efficient  than  the  steady  flow.  In  the  case  of  the 
float  zone,  heat  transfer  from  the  surface  to  the  walls  appears  to  be  more  efficient  with 
detached  oscillating  rolls. 

The  instantaneous  streamline  patterns  in  Fig.  25  support  the  scenario  described. 
Something  novel  is  being  observed:  Rolls  near  the  free  surface  can  extend  toward  the  axis 
and  can  spawn  weak  rolls  within  a  single  closed  streamline.  The  weak  rolls  are  indicated  by 
one  or  two  streamlines  only.  It  should  be  emphasized  again  that  the  streamlines  do  not 
represent  equally  spaced  incremental  values  of  the  mass  flux.  It  is  noteworthy,  that  the 
streamline  patterns  cannot  expose  the  subtleties  of  the  core's  vorticity  field  that  reveal  the 
mechanism  of  oscillation. 

Isotherms  are  given  in  Fig.  26.  The  oscillatory  change  of  the  temperature  field  affects 
the  whole  float  zone,  including  the  axis  r  =  0.  The  temperature  oscillation  in  the  upper  half 
is  clearly  governed  by  convection.  In  contrast,  fluid  motion  near  the  axis  is  virtually  absent, 
and  the  temperature  change  can  only  take  place  by  diffusion.  The  patterns  of  isotherms  next 
to  the  cold  walls  reveal  that  the  change  of  the  temperature  gradient  along  the  walls  is  larger 
than  that  for  a  steady-state  flow  and  oscillating  (Fig.  28).  This  greater  change  confirms  the 
statement  made  in  the  introduction  that  the  solidification  of  molten  silicon  in  an  oscillating 
float  zone  would  take  place  in  a  more  uneven  way  than  that  in  a  float  zone  with  steady 
convection.  The  surface  temperature  is  plotted  in  Fig.  27.  During  the  oscillation  period  the 
maximum  temperature  difference  remains  approximately  3°C  but  its  location  switches  back 
and  forth.  The  periodic  change  of  the  slope  of  the  curves  at  the  corner  is  large,  indicating 
periodic  slamming.  Fig.  27  also  shows  that  the  surface  temperature  is  almost  constant  at  the 
symmetry  point  z'  —  0.5  cm,  that  is,  T,  —  —  2.9 °C.  This  value  corresponds  to 

Ma*  =  943.  A  comparison  of  Nusselt  number  with  radiation  number  reveals  that  the  latter 
is  one  thousand  times  smaller. 

Finally,  computations  were  made  with  the  assumption  of  a  flat  free  surface,  that  is, 
with  h  —  R.  This  assumption  was  introduced  at  t  —  26  s  by  abruptly  changing  the  boundary 
conditions  (7)  through  (10)  and  (13)  to 
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Although  the  actual  free-surface  deformation  is  small  (of  the  order  of  10-4  cm),  the 
combined  or  cumulative  effect  of  “symmetrization”  of  the  free-surface  conditions  due  to  the 
flat-surface  assumption  causes  a  damping  of  the  disturbance  (Fig.  30).  Symmetrization 
means  that  the  terms  in  Eqs.  (7)  through  (10)  and  (13),  which  drop  out  due  to  ht  =  0,  are 
also  terms  which  contribute  to  asymmetric  behavior:  Constant  h  —  R,  Eq.  (25),  and 
vanishing  u  —  0,  Eq.  (22),  obviously  have  a  symmetrization  effect.  Although  less  obvious, 
the  same  is  true  for  the  pressure  at  the  singular  points  R  »  1  cm,  z  =  ±  1/2  cm,  where  the 
flat  free  surface  condition  enforces  p  =  <JufR.  On  a  deformable  free  surface,  no  symmetry 
at  and  near  the  singular  points  occurs,  and  the  pressure  fluctuates.  Also  temperature  and 
temperature  gradient,  Eqs.  (24)  and  (26),  are  affected  by  the  assumption  of  a  flat  surface. 
Near  the  singular  points,  the  term  |  |  in  Eq.  (13)  is  larger  than  the  term  |  Tr  | ,  which  is 

zero  exactly  at  the  singular  points.  On  a  fiat  surface,  vanishing  \hzTt\  near  the  singular 
points  make  the  temperature  distribution  more  symmetric.  This  behavior  influences  the 
vorticity  distribution  along  the  free  surface,  Eq.  (24),  in  such  a  way  that  the  location  where 
“  0  moves  closer  to  the  symmetry  plane  z  =  0  and  with  it  Tmax.  The  oscillation  thus 
diminishes.  It  is  concluded  that  damping  appears  to  be  a  combined  or  cumulative  effect  of 
the  symmetrization  of  all  free-surface  boundary  conditions. 

Damping  of  the  disturbance  means  that  the  oscillatory  mode  of  the  fluid  motion 
changes  with  time  to  a  steady  mode.  This  result  confirms  the  findings  of  Kamotani  et  al. 12 
and  Kazarinoff  and  Wilkowski13  that  a  deformable  free  surface  is  a  necessary  (but  not 
sufficient)  condition  for  the  float  zone  to  become  unstable. 

The  Case  Ma  =  30,225 

In  order  to  narrow  in  on  the  critical  Marangoni  number,  a  third  case  between 
Ma  =  10,400  and  50,050,  that  is,  Ma  =  30,225,  was  chosen  and  computed.  The  results 
show  also  oscillations  and  thus  instability.  This  means,  that  the  critical  Maragoni  number 
must  lie  in  the  interval  10,400  <  Ma  <  30,225.  The  frequency  of  the  oscillations  is 
0.22  Hz,  and  Ma *  at  the  center  is  about  Ma*  =  618.  All  other  features  of  the  oscillating 
motion,  described  for  Ma  *=  50,050,  are  observed  here  too  and  are  qualitatively  the  same. 

CONCLUSIONS 

The  two  main  results  of  this  paper  are  (1)  the  existence  of  a  critical  Mu-number 
between  10,400  and  30,225,  at  which  the  steady  flow  Add  becomes  unstable  and  changes  to 
an  axisymmetric  oscillatory  field.  An  explanation  for  this  oscillation  is  given.  (2)  A 
deformable  free  surface,  even  a  minute  one,  is  necessary  for  the  onset  of  instability.  These 
findings  confirm  those  by  Kazarinoff  and  Wilkowski,13  although  quantitative  differences  exist 
with  regard  to  the  value  of  the  critical  Marangoni  number.  It  is  conjectured  that  the  reason 
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for  damping  is  the  symmetrization  effect  of  the  boundary  conditions  for  a  flat  free  surface. 

Both  the  steady-state  solution  at  Ma  =  10,400  and  the  oscillatory  ones  at 
Ma  =  30,225  and  50,050  consist  of  two  toroidal  rolls  along  the  free  surface  and  two  very 
weak  ones  beneath.  The  essential  fluid  motion  takes  place  in  the  upper  half  of  the  float  zone 
near  the  free  surface.  The  transient  phase  after  the  start  of  the  heating  process  is  practically 
finished  in  less  than  two  seconds  whereas  the  transient  oscillatory  motion  needs  more  than  30 
seconds.  For  Ma  =  50,050  the  frequency  of  the  oscillation  is  0.27  Hz,  for  Ma  —  30,225 
slightly  lower,  that  is,  0.22  Hz- 

The  explanation  for  the  onset  of  instability  and  the  subsequent  maintenance  of  axial 
oscillation  is  sought  in  the  transition  from  attached  rolls  to  detached  rolls  with  increasing 
heat  energy  at  a  critical  Marangoni  number.  Oscillation  is  established  from  an  initial 
disturbance  by  pulsating  rolls  whose  restoring  mechanism  is  the  build-up  and  decay  of  core 
vorticity  supported  by  the  oscillating  location  of  highest  surface  temperature.  The  overshoot 
at  the  symmetry  position  is  due  to  the  large  moment  of  inertia  of  the  detached  rolls  and 
overcomes  diffusion  until  a  quasi-steady  state  with  constant  amplitude  is  reached.  The  flow 
quantity  which  is  best  suited  as  an  indicator  for  instability  and  as  a  descriptor  of  the  pulsating 
flow  process  is  w+/r  and  not  oj+. 

The  free  surface  deforms  only  slightly  of  the  order  of  10-4  cm  in  agreement  with 
Kazarinoff  and  Wilkowski’s  result.13  The  velocity  at  the  free  surface  is  very  high  and  reaches 
values  of  11  cm  /s  for  Ma  ~  50,050. 

The  whole  temperature  field  oscillates  and  causes  uneven  heat  transfer  at  the  walls  that 
is  not  desired  for  silicon  crystal  growth.  It  is  thus  confirmed  that  oscillating  float  zones 
should  be  avoided  in  engineering  applications. 
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Fig.  1.  Sketch  of  a  float  zone  with  streamlines  and  ring  heater  after  Schwabe  et  al.  (1978). 
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Fig  3a.  Mapping  of  the  physical  plane  (r,  z')  onto  the  computational  plane  (£,  i j) 


Fig  3b.  ‘Tour-color”  scheme. 


Fig.  3.  Numerical  grid  layout. 
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0  0.5  1  JO  o  0.5  1.0  cm 
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Fig.  4.  Lines  of  constant  u+/r  for  Ma  —  10,400  at  six  different  times.  The  u^/r  -  contours 
are  ...,  "“3,  “*1,  1,  3,  ... . 
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Fig.  6.  Equivorticity  lines  (lines  of  constant  for  Ma  =  10,400  at  t  —  10.5  s.  The  - 
contours  are  ...»  —3,  — 1,  1,  3 . 
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Fig.  8.  Streamlines  of  the  vorticity  field  in  Fig.  4  at  three  different  times. 
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25 


Fig.  13.  Lines  of  constant  w+/r  for  Ma  =  50,030  at  four  different  times  before  and  after  the 
asymmetric  disturbance.  The  w^/r  -  contours  are  ...»  —18,  —6,  +6,  +18,  ...  . 


Fig.  14.  Lines  of  constant  w^/r  in  the  core  regions  of  Fig.  13.  The  difference  in  the  contour 
lines  is  unity  for  /  *  1.5  s  and  2.0  s,  and  0.5  for  t  =*  2.5  s  and  4.0  s. 


I 


0  0.5  1.0  0  0.5  1.0  cm 

z' 

Fig.  18.  Isotherms  for  Mo  =  50,050  at  the  times  t  —  1.5,  1.72,  2.5,  4.0  s.  The  T  —  Ty  - 
contours  are  0.04,  0.08,  ...,  3.28. 
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Fig.  19.  Surface  temperature  for  Ma  —  50,050  of  the  temperature  field  in  Fig.  18. 
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Lines  of  constant  w^/r  of  a  5/4  cycle  computed  for  Ma  *=  50,050  at  six  different 
times.  The  u^/r  -  contours  are  —18,  -6,  46,  418,  ... 
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Fig.  21,  coni. 
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Fig.  30.  Surface  vorticity  on  a  flat-free  surface  for  Ma  «  50,050  (a)  at  z'  *  0.5  cm  and  (b) 
at  z' —  0.0125  cm  as  a  function  of  time  after  the  change  in  the  surface  conditions  at 
t  -  26  5. 
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